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Abstract. On the base of a Feynman-Kac-type formula involving Poisson stochastic 
processes, recently a Monte Carlo algorithm has been introduced, which describes 
exactly the real- or imaginary-time evolution of many-body lattice quantum systems. 
We extend this algorithm to the exact simulation of time-dependent correlation 
functions. The techniques generally employed in Monte Carlo simulations to control 
fluctuations, namely reconfigurations and importance sampling, are adapted to the 
present algorithm and their validity is rigorously proved. We complete the analysis 
by several examples for the hard-core boson Hubbard model and for the Heisenberg 
model. 
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1. Introduction 

Probabilistic techniques, as the Quantum Monte Carlo (QMC) algorithms, provide 
an essential tool to investigate the properties of many-body systems. Basically, with 
these techniques one evaluates functions of matrices by a random walk in the space 
of the matrix indices [1]. Given the Hamiltonian H of the system, the idea is to 
find an appropriate stochastic representation of the imaginary time evolution operator 
U{t) = exp{—Ht) applied to some initial trial state. By using these methods, one can 
obtain, at least in the absence of a sign problem, the ground-state properties with a 
numerical effort that grows with some power of the size L of the system. On the other 
hand, the exact diagonalization of the Hamiltonian would imply an effort exponentially 
growing with L. 

From a more physical point of view, a probabilistic representation provides a dual 
picture of the quantum systems: on one hand, the traditional description in terms of bra, 
ket and operators, on the other hand, a description in terms of expectations of suitable 
stochastic functionals, which are averages over virtual trajectories of the particles. It 
is this mapping with a (in a sense) classical system that allows us to extract quantum 
information by statistical simulations. 

In recent years, it has been proved that the dynamics of a system of quantum 
particles on a lattice admits an exact probabilistic representation [2, 3, 4]. A suitable 
stochastic functional Mno\ which is defined in terms of a collection of independent 
Poisson processes and diffuses from a Fock state no to a Fock state rit, has the property 
that the expectation value E{A4no^6nt^n), taken with respect to the Poisson processes, 
gives the matrix element of U{t) between the two Fock states no and n. In the 
theory of stochastic processes, this probabilistic representation may be regarded as the 
lattice version of the Feynman-Kac formula. We emphasize that in this method no 
approximation is introduced and no "infinity path integral extrapolation" is requested. 
It will be referred to as the exact probabilistic representation (EPR) of the evolution 
operator U{t). The validity of EPR is not limited to Hamiltonian systems: it can be 
used to express any imaginary- or real-time evolution operator U{t) having any finite 
matrix H as generator. 

In Ref. [5] we used EPR to obtain semi-analytical results in the limit t ^ cxd, in 
which a central limit theorem applies. In this paper, we consider EPR at arbitrary times 
within a Monte Carlo approach (EPRMC). 

Two other well known QMC algorithms, namely the path integrals Monte Carlo 
method (PIMC) and the Green function Monte Carlo method (GFMC), have affinities 
with EPRMC and a comparison is mandatory. 

In PIMC, one evaluates the matrix elements of U{t) by using the Trotter approx- 
imation [6]. The operator U{t) is factorized in the kinetic, exp(— Tt), and interac- 
tion, exp{—Vt), terms so that one gets exp{—Ht) = Y[k=i^w{'~Tt/N) exp{—Vt/N) + 
0([T, y]t^/A^^). This approximation leads to a Feynman-Kac formula, in which, as 
in EPRMC, the trajectories in the Fock space are generated only by the kinetic part. 
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exp{—Tt/N). However, in contrast to EPRMC, there are no stochastic times related 
to Poisson processes. The PIMC simulations are performed by evolving the trajectories 
for TV steps. The drawback is that to obtain results corresponding to t/A^ — 0, in which 
the Trotter approximation becomes exact, one must use extrapolation procedures. For 
any finite value of N, the extrapolation becomes unreliable for values of t sufficiently 
large. This is particularly evident in the case of real times {t it), when the matrix 
elements of U have an oscillating behavior with respect to t. In contrast, no small step 
approximation is introduced in EPRMC and no extrapolation is requested. 

Now, let us consider GFMC The method consists in repeated statistical 
applications of the operator (1 — Ht/N) to an initial state. For oo, one 

reproduces the operator U{t), whereas an approximation affected by a relative error 
e{N) is obtained for any finite A^. It is plausible that the sampling directly the operator 
U{t) instead of (1 — Ht/N)^ leads to a higher efficiency. In the Appendix, we show 
that the relative efficiency between EPRMC and GFMC in filtering the ground state is 
Ei/[2E^^\Ei-Eo)e], where Eo and Ei are the energies of the ground- and first-excited 
states of the considered system and E^l^^ is the ground-state energy of the associated 
non interacting system. Since the gap {Ei — Eq) decreases as the size L of the system is 
increased, compared to GFMC, EPRMC offers a more powerful method to investigate 
the ground-state properties of large lattice systems. 

Actually, the GFMC scheme mentioned above is rather crude. Trivedi and Ceperley 
[7] introduced Poisson processes as a tool to obtain a more efficient GFMC method 
when the transition probabilities, proportional to the matrix elements of Ht/N, vanish 
for A^ ^ oo. We will refer to this improved GFMC as GFMCP. In Ref. [4] it has been 
proved that in the limit A^ ^ oo GFMCP becomes equivalent to EPRMC. However, as 
explained in the Appendix, for a finite value of A^ the relative efficiency of EPRMC with 
respect to GFMCP is {Eq/ e'^'^Y /2e, i.e. it is proportional to the accuracy required 
in the approximated GFMCP. 

Controlling the large fluctuations is one of the most important issues of any Monte 
Carlo method. This is evident in GFMC where an iterated statistical application of 
the operator (1 — Ht/N) is performed. Roughly speaking, after k iterations one has 
fluctuations that grow like A*', A being the statistical error associated to a single step. 
To solve the problem of large fluctuations, besides the development of the importance 
sampling method [6], remarkable progress has been made with the reconflguration 
technique introduced by Hetherington [8] and subsequently improved by Sorella [9] who 
proposed a scheme without bias (see also Ref. [10]). 

In this paper, after introducing some relevant physical models (Section 2) 
and recalling EPR (Section 3), we extend EPR to the study of exact time- 
dependent correlation functions (Section 4). In the core Section 5, we discuss the 
EPRMC algorithm flrst with a pure sampling and then adding fluctuation control 
by reconflgurations and importance sampling. A detailed proof of the validity of the 
reconflguration method is given in Section 6. Results of numerical simulations for the 
hard-core boson Hubbard model and for the Heisenberg model are discussed in Section 
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7. Conclusions are drawn in Section 8. 
2. Models 

The Hamiltonian models of interest have the following general structure (we shall always 
take h = 1) 

H = T + V, (1) 

where V is the potential energy operator and T the kinetic energy operator, which on 
a lattice assumes the form 

^ = - E E (^Is^ + 4^-) • (2) 

Here A C 2''^ is a finite d-dimensional lattice with | A| ordered sites and Cjo- the commuting 
or anticommuting destruction operators at site i and spin index a with the property 
cf^ = (fermion or hard-core boson systems). The system is described in terms of 
Fock states labeled by the configuration n = (nij^, riij^, . . . , ?t,|a|i-, ''^|A|i), ^-e. the set of 
lattice occupation numbers taking the values or 1. The total number of particles is 
N„ = XliGA for cr =t J,. We shall use the mod 2 addition n (B n' = [n + n') mod 2. 

The analysis we develop in the following is valid for an arbitrary functional form 
of the potential V . However, numerical examples will be limited to the well known 
Hubbard potential [11] 

V = Y^lic\^Ci^ c\^Ci^. (3) 

ieA 

We emphasize that, independently of its form, V is diagonal in the Fock space, whereas 
T is off diagonal. 

In this paper we will consider only hard-core boson systems. We recall that, even 
if fermion systems, like the Hubbard model, are more attractive, hard-core bosons have 
not a purely academic interest. Besides the description of boson particles with a hard- 
core interaction, they can be mapped onto systems of half integer spin [1, 7, 12]. As an 
example, we consider the S = \ Heisenberg quantum antiferromagnetic model 

H=JY,S^-S, = i^ Y.{StSj + SrSp + jY^StS',. (4) 

where J > 0, (i, j) indicates that the sites i and j are nearest neighbors, and Si and Sj 
are the spin operators. It is convenient to view the left and right factors in Si ■ Sj as the 
spin operators of two sublattices A and B, respectively. The mapping is then established 
as follows. The operators 5*^^ and S~ commute on different sites and are thus identified 
with boson operators via b] = and bj = S~. Furthermore as S- = S^S~ — |, 
one has = rij — i, where Ui = b\bi is the number operator. For a half spin system 
SfSf = S~S~ = 0, which implies = 0. Therefore, the bosons have a hard core and 
a site can be occupied by at most one particle. In order to a have negative sign in the 
kinetic energy term, a further transformation is necessary. The spins on the sublattice B 
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are rotated as Sj ^^j^ ~^ ^ ^^"^ ~^ ^j- Since this transformation is unitary, 
the commutation relations are left unchanged. The hard-core boson Hamiltonian then 
reads 

where Ej\f = —JZ\A\/8, Z being the number of nearest neighbors for the given lattice, 
e.g. Z = 4 for a square lattice in two dimensions. 



3. Probabilistic representation 

We are interested in evaluating matrix elements of the form (n|e~'^*|no) or (n|e~*^*|no) 
between two Fock states Uq and n with t G M. As usual, we will speak of imaginary 
times in the former case and of real times in the latter one. 

Let r be the set of links, i.e. the pairs (i, j) with I < i < j < such that rjij ^ 0. 
For simplicity, we will assume rjij = r] Hi and j are first neighbors and r]ij = otherwise. 
For a d-dimensional lattice the number of links per spin component is |r| = d\A\. Let 
us introduce 

Xij^in) = {n® li„ © IjalcLs'. + c]^cjn), (6) 
V{n) = {n\H\n), (7) 

where lia = (0, . . . , 0, lia, 0, . . . , 0). Note that the values assumed by Xija are or 
1 (Ajjcr = 0, ±1 is possible in the case of fermion systems not considered here). We 
will call the link (ija) active if Ajjo- 7^ 0. Let {N^j^-}, ihj) G F, be a family of 2|F| 
independent left continuous Poisson processes with jump rate p = r] ii Xij^ ^ and 
otherwise [13]. Let us now define the stochastic dynamics on the lattice. At each jump 
of the process N\^^ a particle with spin a moves from site i to site j or vice versa. Let 
us indicate with A[n) the number of active links in the configuration n 

(ij)Gr<7=Ti 

The total number of jumps at time ^ is A^t = j)er X]cr=tj, ^ija- By ordering the jumps 
according to the times Sk, k = 1, . . . , Nf, at which they take place in the interval [0, t), 
we define a trajectory as the Markov chain ni,n2, . . . ,nN^ generated from the initial 
configuration no by the stochastic dynamics described above. Let us call Ai, A2, . . . , Aat^, 
Vi,V2 ... ,VNt and Al,A2...,AN^ the values of the matrix elements (6), (7) and (8) 
occurring along the trajectory, respectively. For simplicity, we will indicate the last 
configuration reached after Nf jumps as rif = n^^. We will also use the symbols 
Ao = A{no), Vo = V{no) and sq = 0. 

According to Ref. [4], the following representation holds 

(n|e"^*|no) = E(5.,.,A^Sf), (9) 
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where the expectation E (■) has to be taken with respect to the Poisson processes N^^^ 
and the stochastic functional TWSq*^ is defined by 

The subscript no in Ain'o '' specifies the initial state. For real times an analogous 
representation holds 

(n|e-^*|no) = E(Vn.A^lSf), (11) 

where 

J\yl[0,it) ^ -Nt ^Jo[A(ns)v-iV(n.s)]ds ^ ^-^2) 

In the following, we will consider the case of imaginary times. Except when explicitly 
said, all the formulas are trivially extended to the case of real times. 

Any ground-state quantity can be obtained from the matrix elements (9) by a 
proper manipulation and taking the limit t ^ oo. For instance, the ground-state energy 
is given by 



-Er^^t{ 


n 


e ■ 


Ht\ 


l^o) 




e" 




l^o) 



= lim ' °^ = T > . (13) 



It is easy to see [4] that —dtE{M.no ) = E.{Mna H{nt)), where H{nt) is given by 

Hint) = J2{^'\H\nt) = -[A{nt)v - V{nt)]. (14) 

n' 

Equation (14) is the local energy of the last visited configuration n^. Therefore, Eq. (13) 
becomes 



E(7^(r^,)^ '°'*^ 



'■no 



These identities are valid if the initial configuration tiq is such that (-EqI^o) 7^ 0. For a 
finite t, this scheme allows a good estimate of Eq ii t ^ {Ei — Eq)"^, where Ei is the 
first-excited state of H. This implies that t must be increased by increasing the lattice 
size |A|. 

4. Correlation functions 

Let us consider a generic operator O. By using twice the Fock representation of the 
identity operator and twice Eq. (9) with functionals Ain'o^ and A^'^°'* \ respectively 
defined by two sets of independent Poisson processes {N-j^} and {A"/*'^}, we have 

(n|e-^*'Oe-^*|no) = J2J2^n\e-'''' \n'^){n',\0\n"){n"\e-'''\no) 

n'a n" 

= Y: E (5.;,,. M'Sf^ {n',\0\n,) A^L?) , (16) 
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where n'^, is the configuration reached at time t' starting from tIq. From this expression, 

we get 

J](n|e-^*'Oe-^*|no) = E ^ (-M'^ ^ri|0|n,) M^^f) . (17) 

n n 

The ground-state quantum expectation of the operator O, assuming (i^ol-Eo) = 1; iS) 
therefore, 



{Eo\0\Eo) = hm 



t, t'^oo ^ ( _A/f[°'*+*') 



I no 

We can consider two basic cases for the operator O. 
4-1. Diagonal operators 

In this case, {n'\0\n) = 6n'.nO{n) and Eq. (18) becomes 



E I^mT'^ Oint) M 



no 



(Eo|0|i?o)= hm ^ ^. (19) 

Note that E{M'!^P Mtf) = E(A^[2;/+*'^), so that {Eo\0\Eo) = 1 if O is the 
identity operator, whereas for a single reahzation of the stochastic functionals we have 

mT'^ A^iSr'^ 

4-2. Off-diagonal operators 

In this case, O is typically given in terms of elementary operators Oija connecting 
two different Fock states like (n'|Ojjo-|n) = Ojjcr(n)(5„/ ,fji<T«i<T, where n^"'^^'' is the 
configuration obtained from n exchanging Uia- with Uj^. Therefore, one has 

(i^^o|0.,.|i^o) = lim ^ . ^. (20) 



Similar expressions hold for other off-diagonal operators connecting two generic Fock 
states. 

5. EPRMC algorithm 

5.1. Pure sampling 

Equations (9) and (10) lend themselves to a statistical evaluation of the matrix elements 
(n'|e~'^*|n) via a random sampling of jump times and trajectories. As explained in 
Ref. [4], the practical algorithm works as follows. We start by determining the active 
links in the initial configuration no assigned at time and make an extraction with 
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uniform distribution to decide which of them jumps first, say the hnk We then 

extract the jump time Si according to the conditional probabihty density 

PAois) = Aor]exp{-Aor]s), (21) 

where Aq is the number of active hnks before the first jump takes place. The contribution 
to Ain'o^ at the time of the first jump is, therefore, 

^{Aov-Vo)s^Q^^ _ ^ ^(Aov-Vo)tQ^^^ _ (22) 

According to Eq. (10), the contribution of a given trajectory is then obtained by 
multiplying the factors corresponding to the different jumps determined in an analogous 
way until the last jump takes place later than t, i.e., 

_/\/([0,t) = ■J~|' g(^fc-i»?-Vfc_i)(sfc-Sfe_i) I ^{ANtV-VNt){t-SNt) (^23) 

if Nt > 0, and M^^f = e^^o^'-^o)* if Nt = 0. 

Let us consider independent trajectories obtained as described above and let 
Mn'o^^^^ be the functional value (10) calculated along the 2-th trajectory. From the law 
of large numbers we have 

1 ^ 

E(A<£*^)=JL- (24) 

i=l 

5.2. Reconfigurations 

Equation (10) represents a product of Nt random factors and, since Nt grows with t, 
the fluctuations of the functional tMISq*'' grow exponentially with t. This implies that 
the number of trajectories needed to have good statistical averages grows exponentially 
with t. A similar problem has been successfully tackled some years ago in the framework 
of GFMC by the reconfiguration technique [8, 9]. This technique can be adapted also to 
the present probabilistic representation. In fact, for boson systems at imaginary times 
the stochastic functional TVln,*'' is always positive and can be thought as a weight. Let 
us divide the time interval [0,t) in R subintervals of the same length At = t/R. Let us 
label the times corresponding to the end-points of these intervals as 

tr = rAt, r = 0, (25) 

and let n^,. be the configuration reached at the time ^,.-1-0^ trough the dynamics described 
in Section 3 (we recall that the Poisson processes are left continuous defined). The 
following obvious identity follows from Eq. (10) 

R 

A^l?;*) = n-^fe;_f^ (26) 

r=l 

which implies 

E(A<£*^)=E^n-^t_f)V (27) 
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The functional Ai^nt^^f'^^ will be referred to as local weight. 

Essentially, the idea of the reconfiguration technique is the following. Instead 
of extracting independent trajectories, one carries on an ensemble of M trajectories 
simultaneously in order to perform dynamically, at the times tr, a suitable replica 
of those with large weights, eliminating at the same time the others. This 
replication/elimination of trajectories, also referred to as reconfiguration, has to be done 
in such a way that the number of trajectories M remains constant. At the end, one can 
substitute the average ofRti-^-v-f^with the average 

rit^ rit^ we indicate the reconfiguration action at the time t^? and with (A^-''"^'*'''') 
the uniform "average" of the local weights over the M reconfigured trajectories (we use 
quotation marks since this quantity is itself a random variable). Hence, the remarkable 
advantage of using the reconfigurations is that, if the functional H^Li -^ntjT^f'^'* has 
variance A^*, the variance of ri^Li(-^li'^~^'*'^'') "^i^ be roughly (A/-\/M)^*, where A is 
the variance of the local weights and R* < R is the number of subintervals in which the 
local weights become approximately independent. 



5.2.1. Reconfiguration algorithm, here, we describe in detail the reconfiguration 
algorithm at imaginary times postponing the relative proof to the next Section. We 
will indicate with rif \ r = 0, . . . , i?, and _A/i[*^-i'*'")(*)^ y = 1, . . . , i?, the configurations 

and the local weights of the z-th trajectory and define the corresponding M-component 
vectors as n^^ and M}nt^]_'^'^\ respectively. We shall use also the operator symbols V 
and TZ: V applied to the configurations rit^ gives the configurations rit^^^ according to 
the dynamics defined in Section 3 along the time interval [t^, ^r+i), whereas TZ gives the 
reconfigured configurations n^^ = TZrif^.. 

First step. Define ritg = rit^ with ■n!^^ = uq for i = 1, . . . , M. At the initial 
time to = 0, all the M trajectories starting from the initial configuration no follow the 
dynamics V and reach the configurations ^ = T>ntQ . Correspondingly, evaluate the M 
local weights along the time interval [0,ti), and compute their average 

^0 

M 

(A^Lo.*0^^^y_^[o,t.)(0_ (28) 
1=1 

Second step. Perform the reconfiguration rit^ rit^ = TZrit^. The new 
configurations are obtained by drawing out them randomly from the old ones, n^j, 
according to the probabilities 

_y^[0,ii)(i) 

' ti - , .[o,ti){0- ^ ^ 

2^1=1 ■'^^ntg 

The new configurations n^^ are used as starting configurations of the M trajectories 
for the time interval [^1,^2) and, through the dynamics V, are mapped into 'Drit^. 
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Correspondingly, evaluate the M local weights A^~^'*^^ and compute their average 



M 

-FE-^^if""- (30) 

Third step. Perform the reconfiguration Vrit^ nt.-, = IZVnt-^ by drawing out the 
new configurations randomly from the old ones according to the probabilities 

The new configurations rit^ are used as starting configurations in the time interval [^2, ^s)- 
Evaluate the local weights tW-^'*^'' and compute their average 



(32) 



M 
1=1 



By iterating this procedure for R steps, we arrive to the final configurations 



Vnt^_^ = V(nV)^^^nt^, with R computed averages (.mI^''*''^), r = 1,...,R. As 
we will prove later, the following identity holds 

R 

l[tr-l,tr)\ 



E (-MlS^) = E n(-^t:f ^) ' (33) 



where E indicates the expectation in which the configurations n^^ are obtained by the 
reconfiguration procedure described above. Explicitly, Eq. (33) implies that to evaluate 
the expectation E(7WlSo*''), instead of Eq. (24), we can use 

E(A<£f)= hm f[{Mtr':^), (34) 

r=l 

or, more generally, simulating independent samples each one composed by M 
reconfigured trajectories, 

, N R 

E(A<£") = ,to^^En(-^£;';''')'"- ps) 

p=l r=l 

The label (p) in Eq. (35) means p-th sample. Note that for M = 1 we recover Eq. (24). 

All what we said about the functional Ain'o'' can be repeated for the functional 
A4no^Sn^n^- In this case, Eq. (33) becomes 

\r=l 1=1 '-f^-i / 

Equation (36) allows to calculate the numerator of Eq. (15) as 

E {M^^fnin,)) = E n(-^tf ^)]^ "*^^'^^((^^)^'^) • (37) 

\r=l 1=1 J 
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5.2.2. Correlation functions. Let us now consider the reconfiguration procedure for the 
functionals introduced in Eqs. (19) and (20) to obtain the correlation functions. In this 
case, we perform R steps in the first interval [0, t) and R' steps in the second interval 
[0,t'). All the quantities relative to the second interval [0,t') will be be indicated with 
a prime. In the pure sampling, the initial configurations of the second interval [0, t') are 
equal to the final ones of the first interval [0,t): n[,^ = ritp,. For diagonal operators, we 
have 

r=l r'=l 

M \ 

><7^E-^5r"*'''o((^^'^-*-)^') ' (38) 

1 = 1 J 

where now the configurations TZ^'Vutj^^ are obtained by updating the intermediate 
configurations at time tn, namely Vnt^_-^, R' times according to the successive R' steps. 
For off-diagonal operators, we have 

r=l r'=l ''■'-1 

M \ 

X ]^ , (39) 

1=1 V-i / 

where , r' = 1, . . . , i?', are the configurations obtained after r' steps starting from 

r' — 1 



the intermediate configurations iT>rLt^_^y^^^" in which the occupations of sites i and j 
with spin a have been exchanged, i.e. n^f = {TZ'Dy'~^{T>nt^_^y'^'^^'^ . 

r' — 1 

5.2.3. Real times. A reconfiguration procedure can be performed also at real times. In 
this case, the stochastic functional A^l^o**'' is complex and we separate the contributions 
from the R time intervals in their moduli and arguments, i.e. 

■ I [*r— 1 i^r ) 



where 



^irZt^ = liN,,. - N,^_,) - V{nMs. (42) 

The moduli can be used as local weights for the reconfiguration operator TZ. All the 
steps described in Section (5.2.1) remain unchanged except for the last factor, which 
takes into account the R phase factors reconstructing the original stochastic functional. 
The final result result is 

E (-MlSf ) = E (U(\^'iT'\) ^ E \M^'"'''''\ e ) . (43) 



,r=l 1=1 



^R-l 
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5.3. Importance sampling 

Although the reconfiguration method controls the growth of the fiuctuations of TW'^'*) 
along the trajectories, since the dimension of the Fock space grows exponentially with 
the lattice size, an extraction of the jumping links by importance sampling also may 
be mandatory to reduce the statistical errors of the local weights [6]. If some a priori 
approximation \g) of the ground state is known, which has the property {n\g) G M \ 
for any Fock state |n), then instead to sample directly the operator exp{—Ht), it 
can be notably advantageous to sample the isospectral operator exp{—Hgt), where 
{n'\Hg\n) = {n'\g) {n'\H\n) {n\g)~^ . 

As explained in Refs. [4] and [14], if \g) is a guiding function in the sense specified 
above, the generalization of the present algorithm to the case with importance sampling 
consists in replacing the number of active links, A{n) = J2{ij)£r J2a=u l-^yo-l^)!; 
the previous formulas with the quantity 

{n ® ® lj^\g) 



(44) 



{n\g) 

Correspondingly, the probability density for the jump times becomes 

Pa,{s) = Agr]exp{-Agris), (45) 

and the extraction of the jumping link {i,j, cr) among the active ones must be performed 
according to the probabilities \{n © Ijo- © ljcr\g){n\g)~^\/Ag{n). Finally, the stochastic 
functional (10) is modified as 

The advantage of using importance sampling becomes clear considering the local 
energy associated to Hg 

Hgint) = y^{n'\g){n'\H\nt){nt\g)-' = -[Ag{nt)v " ^K)]- (47) 



In fact, in the limit \g) \Eo) one has Hg^rit) Eq and accordingly Ai^g'^ — >■ 
exp(— E'ot) so that the fiuctuations vanish. 

For any choice of the guiding function \g), the modified stochastic functional 
(46) provides unbiased representations of the ground-state energy Eq as well as of the 
expectation of a generic operator O in the ground state \Eo) of H. In fact, Eq. (15) now 
reads 

lin, . (48) 

E(XK) 

where Egg is the ground-state energy of Hg, which, however, is an operator isospectral 
to H. On the other hand, Eq. (18) written in terms of a gf-modified operator Og becomes 



lim ^ / ^owM^ '- = {Eo,\Og\Eoo)- (49) 
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By using (njiJog) = (^|5')(*^|-E'o) and {EQg\n) = {n\g)~^{EQ\n), it is simple to see that 
{EQg\Og\Eog) = {Eo\0\Eo) H we choose Og as the operator defined by {n'\Og\n) = 
{n'\g) {n'\0\n) {n\g)~'^ . Note that, in the case of diagonal operators, Og = O. 

Importance sampling may be useful also for a different purpose, namely the 
determination of the transition amplitudes {g\e~''^^\no) between two chosen states |no) 
and \g). This is particularly interesting at real times and we illustrate the idea in this 
case. If \g) is a generic state with the property (n|(7) G M \ so that the isospectral 
Hamiltonian Hg is well defined, we have 

5^(n|e-^«*|no) = (^?|no)((7|e-^*|no). (50) 

n 

Since the expectation of the stochastic functional Ain'o^'^ with the modified rules (44) and 
(45) provides an exact representation of the l.h.s. of Eq. (50), we obtain the transition 
amplitudes (5'|e~*^*|no) up to the constant {g\no). 



6. Proof of the reconfiguration algorithm 

In this Section, we prove Eqs. (33-39). Let us consider an ensemble of M simultaneous 
trajectories obtained by the dynamics described in Section 3 starting from the 
initial configuration no. Let PR{Mn't^^\ TWll'f , . . . , -A/lltf^_'f ^ ; rito, rh^, nta) be 
the probability density to have a realization in which the M trajectories have local 
weights J^^nto^\ J^^nt'^'^\ . . . , J^nf^^f and configurations n^g, rifj, , . . . , n^^ at the times 

to, ti, . . . , t_R, respectively. For simplicity, here we shall often use 

and rir for rit^. Since the M trajectories are independent, if we take n^^ = no for 
/ = 1, . . . , M, we have 

E (A^kr^r.,..) = E (f[Mr^,5^,^)j = E (^f:flM^:l^S^^^i^ • (51) 

\r=l / \ 1=1 r=l J 

Consider, then, the following identity 



1=1 r=l \l=l / \l=l / \l=l 



2 



X 



M \ 

where the quantities po,pi, . . . ,pr-i are defined recursively by 

(i) 1 

i^) _ M^J^, • (53) 

~ ll^UM^l.P^l^ r-i,...,i< i 
Equations (51) and (52) lead to 

E (^Mr^l5r.,r.A = E (^{Mr-l) ^{M R-l5r,,r^a) v\ , (54) 
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where the weighted "averages", {J^r)w and are defined as the weighted 

sums {Mr)w = Y^fii-^r^Pr^ ^nd (A^ij_i5n,rifl)u, = T^fLi -^R-^i^ ^^^wPr-i, respectively. 
Up to now the quantities pr have been thought as stochastic variables. Actually, 

since the components are positive and normalized to 1, we can interpret them 
as probabilities to modify the original probability density Pr. We introduce a 
new probabihty density Pr that, besides taking into account the dynamics V^, 
includes the probabilities Pr, for r = 0, ...,-R — 1. In this case, if we indicate 
with no, ni, . . . , n/}_i, Dn^.i the configurations extracted according to the probability 
density Pr, Eq. (54) transforms into 

(R \ _ /-^"^ 1 \ 

n-^-i^^J = E [U{Mtlf^)-J2M^^^^^^^ , (55) 

where E(-) means expectation with respect to Pr and the weighted "averages" {Air)w 
have been substituted by uniform "averages" over the new configurations, {-M-^'^f'^) = 

"•r-l 

Equations (51) and (55) reproduce Eq. (36). To conclude the proof, we still have 
to show that the algorithm described in Section 5.2.1 coincides with sampling the 
configurations no, ni, . . . , n/j._i, X'n/j„i according to the probability density Pr. For 
M trajectories with local weights A^|,*li, let us define the following probabilities 

Due to the recursiveness of Eq. (53), for r > 1, we have 



p 



(*)-an^r'^ (57) 



r 



where Cr is a normalization constant independent of the trajectory index (i). This allows 
to realize the transformation Pr Pr recursively for r = 1, . . . , i?. At the first step 
r = 1, since po is uniform we do not have to reconfigure and uq = tiq. The density Pi 
will be then sampled through the vectors n-o and Duq. Suppose now to have sampled the 
density Pr through the vectors no, ni, . . . , n^-i, Vrir^i. To sample the density P^+i we 
must change the arrival vector of configurations Vrir^i into a new vector n^ according 
to the probabilities Vr, with components 

With a further dynamic step we get Vrir. The distribution is sampled by iterating 
this procedure R times. This is exactly the procedure explained in Section 5.2.1 and 
the reconfiguration algorithm is proved. 

Equation (33) follows easily by summing Eq. (36) over n. Finally, Eq. (37) can be 
obtained multiplying Ain'o^dn^m by Ti.{n) and then summing the product over n. 
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Let us now consider the functional Ai'n/ Sn^nt-^n'o'^ ■ analogy to the previous 
case, we easily arrive to 

(R R' \ f ^ \ 

n-^r-l^",-H n ^r'^l = E \{{Mr-l)^ n {K'^l)^{^R'^l^r.,r^^)^ , (59) 
■r=l r'=l / \r=l r'=l / 

where, recalling that n[,^ = n^^, the weighted "averages" are given in terms of 
probabilities pr defined recursively as in Eq. (53) with r = l,...,i?+_R' — 1. Let 
Pr+r' and Pr+ri be the obvious extensions of the distributions Pr and Pr previously 
considered. As before, by using Eqs. (56) and (57), for r = 1, . . . , i? + i?' — 1, we can 
realize the transformation Pr+ri Pr+r' recursively: along the interval [0, t) we sample 
Pi, P2, • • • , Pr, whereas along [0, t') we sample Pr+i, Pr+2, ■ ■ ■ , Pr+r', obtaining the 
configurations no, ni, . . . , Ur, Uq, n[, . . . , n'^i_^,'Dn'^,_^. Therefore, Eq. (59) transforms 
into 

R R' \ / R R'~l 



\r=l r'=l 
1 ^ \ 

i^E-«!".'r"""""^.,««'«.-,)<» . (60) 



\r=l r'=l / \r=l r'=l 

M 

i ■ 

X- 



1=1 ^ ^ 



which yields Eq. (38) after multiplying Aino^6n,nt by 0{n) and then summing 
over n. Note that in the r.h.s. of Eq. (60) there appears TZ^'Vur^i and not 
VriR^i. Indeed, according to Eq. (52), in the last weighted average {Ad'R,^i5nnR)w = 

there appear the probabilities p'^^^/_i associated to the last 

time interval. 

In general, in the reconfiguration procedure a weighted "average" 

(-^llrV*''V(rio, ni, . . . , n^_i, ur))^ (61) 
will be substituted by the uniform "average" 

M 

i- J^M^J^'-^^^Hn^-'nofK {n^-'n,f\ . . . , {nnR_,f\ {VuR^fyy (62) 
1=1 

Equation (39) can be obtained in the same way as Eq. (38). Finally, in the case of real 
times, Eq. (43) is immediately obtained by using for the local weights the quantities 
l-^jT.t^Ii'**'^'*! and for the function /(■) of Eq. (62) the product of the phase factors 

/=J]e'*"v-i . (63) 



7. Numerical results 



In this Section, we present some numerical applications of the algorithm described 
above. In principle, the reconfiguration scheme can be applied for any positive integer 
R. However, we have observed optimal reconfiguration for R ~ {A)pt, where (A) is 
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the average number of active links. This is what one expects as, in this case, the 
reconfiguration is repeated in the average at each jump, i.e. as frequently as the 
stochastic dynamics dictates (see also Ref. [15]). In the simulations reported below, 
therefore, we always work with this approximately optimal number of reconfigurations. 

The count of the active links and of the potential of a given configuration, quantities 
to be determined at each jump, is a core point of the algorithm. Starting form a first 
count based on a systematic inspection of the initial lattice configuration, we have 
implemented a local updating of these quantities. In fact, when a jump occurs, the new 
Hubbard potential and the new number of active links are determined by the change 
of the lone occupations and of the lone links involved in the jump. The computational 
cost of a local update, which takes into account only these relevant sites and links, is 
independent of the lattice size. Also the reconfiguration procedure has been optimized 
by defining a non-negative integer, the replication multiplicity where (i) is the 
trajectory index and ^fti l^r^ = Configurations for which jjif' = are substituted 
by those with /^l*^ > 1, whereas no operation is performed for the trajectories with 
/ir ^ = 1, which are the largest fraction of the whole set of M trajectories. The efficiency 
of the resulting code can be figured out by the following example. With an ordinary 
personal computer and without using importance sampling, we are able to simulate 
lattices with 40 x 40 sites with 800 hard-core bosons obtaining the ground-state energy 
up to a relative error of the order of 1% with 290 minutes of cpu time. A detailed 
comparison of the efficiency of our EPRMC code with those implementing other Monte 
Carlo methods is beyond the purposes of present work. In the Appendix, we discuss the 
relative efficiency between EPRMC and GFMC or GFMCP. 

In Figs. 1-5 we compare several quantities evaluated by the EPRMC algorithm with 
the corresponding exact results obtained by numerical diagonalization of the associated 
Hamiltonian. The system considered is a hard-core boson Hubbard model of small size, 
namely a 2 x 3 lattice at half filling. The general purpose of these figures is to show the 
unbiased statistical convergence of the Monte Carlo data towards the exact values. No 
importance sampling is used in these first examples. 

In Fig. 1 we show the expectation E{Aino^) as a function of the imaginary time 
t. The agreement with the corresponding quantum matrix element ^„(?^|e~^*|no) is 
excellent. The reconfiguration procedure is able to control completely the fiuctuations 
growing with t so that the error bars, negligible on the used scale, do not increase by 
increasing the time. Similar results are obtained for different initial configurations no. 

In Fig. 2 we show the expectation E{A4no*^) as a function of the real time t. 
Also in this case there is an exact statistical convergence towards the quantum matrix 
element ^„(?T'|e~*^*|no). However, in this case the reconfiguration procedure is able 
to control only a part of the fiuctuations, namely those related to the modulus of the 
functional Mno^\ The fiuctuations associated to the corresponding phase factor make 
the convergence harder and harder for large times. 

In Fig. 3 we show the behavior of the local energy E{Aino^H{nt)) /E{Aino^) as a 
function of the imaginary time t. According to Eq. (15), the local energy converges 



Exact Monte Carlo time dynamics in many-body lattice quantum systems 17 

,80 




Figure 1. Expectation of the functional Mn'o ^ versus the imaginary time t for 
a hard-core boson Hubbard system in a 2 x 3 lattice at half filling with rj = 1, 
7 = 4, and periodic boundary conditions. The initial configuration is no = 
(1,1,1,0,0,0,1,1,0,1,0,0). The Monte Carlo simulation (dots with error bars) was 
done with M ~ 2^^ trajectories, N = 2^ samples, and R = 300 reconfigurations. 
Error bars correspond to one standard deviation evaluated from the N samples. The 
dashed line is the exact result from numerical diagonalization of the corresponding 
Hamiltonian. In the inset we show the small time behavior. 



to the ground-state energy of the system, Eq, for large times. In fact, after an initial 
transient inversely proportional to the gap Ei — Eo, the ratio E{Aino^H{nt))/E{Aino^), 
estimated with a finite number of trajectories M, fluctuates around an average value 
that is close but different from Eq (see inset of Fig. 3). However, if we increase M, as 
shown in Fig. 4, the statistical accuracy increases and we obtain an unbiased convergence 
toward Eq. 

As an example of correlation functions, we studied the spin-spin structure factor 
S{q., qy) = m E e'''^^^^-^^^+"^y^y^-y^\E,\S,S,\Eo), (64) 

where Si = c^^c^ — cJjCj| and Xi and yi are the coordinates of the i-ih. lattice point. Note 
that the operators SiSj are diagonal in the Fock space and can be evaluated by using 
Eq. (38). In Fig. 5 we show S{'k^tt) evaluated for different values of the interaction 
strength 7. In agreement with the exact results from numerical diagonalization, S{tt, it) 
has a sharp transition between the 7 — > and 7 — > 00 asymptotic values. This transition 
is expected to take place when the average kinetic and potential energies are of the same 
order, i.e., for r]{A) ~ 'j{N-^i), where (iV^j.) is the average number of doubly occupied 
sites. For the system considered in Fig. 5, we have (A) ~ 15 and (A^'^j,) — 1.5 so that 
the transition is expected at 7/77 ~ 10. This is in agreement with the numerical results. 
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Figure 2. Expectation of the real and imaginary parts of the functional 
versus the real time t for the same system of Fig. 1. The Monte Carlo simulation (dots 
with error bars) was done with M = 2^° trajectories, N ~ 2"^ samples, and R ~ 15 
reconfigurations. Error bars correspond to one standard deviation evaluated from the 
N samples. The dashed (real part) and dot-dashed (imaginary part) lines are the exact 
results from numerical diagonalization. 




-10.5 



Figure 3. Local energy E{Aina ^'H{nt))/E{Mni;^^) versus the imaginary time t for the 
same system of Fig. 1. The Monte Carlo simulation (solid line) was done with M = 2^'', 
= 1, and R ~ 300. The straight dashed line is the exact energy Eq = —10.233803 
obtained by diagonalization. In the inset we evidence the difference between Eq and 
the time average of E{Mno^Ti.{nt))/E{Mno'^) computed over the interval 10 < t < 20 
(opaque region baseline). 
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Figure 4. Relative error between tlie local energy E(xl?o*^H(nt))/E(M[?^*^) and the 
exact energy Eq versus the number M of reconfigured trajectories for the same system 
of Fig. i with = 2^, < = 5, and R = 75. Error bars correspond to two standard 
deviations evaluated from the N samples. 




Figure 5. Spin-spin structure factor S{qx,qy) at Qx = Qy = tt versus the interaction 
strength 7 for the same system of Fig. 1. The dashed line is the exact result from 
numerical diagonalization of the Hamiltonian whereas the dots with error bars are 
from a Monte Carlo simulation with M = 2^^ {M = 2^^ for 7 > 10), N = 2'^, 
t = t' = 5, and R — 4:5. Error bars correspond to one standard deviation evaluated 
from the N samples. 
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Figure 6. Local energy per site [E(7W[?o*^H(nt))/E(7W[?Q*^)]/|A| versus the imaginary 
time t for two different-size hard-core boson Hubbard systems at half fiUing with rj ~ I, 
7 = 4, and periodic boundary conditions. The Monte Carlo simulations (solid lines) 
were done with M = 2^^, R = 2^^, and TV = 1. The straight solid lines are the time 
averages of the simulation results computed over the interval 10 < i < 40 whereas the 
straight dashed lines indicate the relative standard deviations. The simulations took 
79 (20 X 20 lattice) and 290 (40 x 40 lattice) minutes on a computer with a 2.40GHz 
Intel Xeon CPU. 



In Fig. 6 we report simulations performed for hard-core boson Hubbard systems of 



large size. In particular, we show the local energy per site [E{Aino^H{nt))/E{Aino')]/\A\ 
as a function of the imaginary time t for two lattices at half filling having size 20 x 20 
and 40 X 40. Note that the standard deviations of the fluctuations around the long-time 
averaged value of [E{M.no^H{nt))/E{Aino^)]/\A\ provide an estimated relative error for 
Eq/\A\ of the order of 1%. This result is obtained with a moderate computational effort. 
In Fig. 6 it is also evident an asymmetry of the fluctuations of the local energy around 
its mean value. This behavior is due to the reconfiguration procedure that ensures the 
invariance of the first statistical moment of A^'^'*^ (or of related quantities) only. 

We have performed simulations also for the Heisenberg model (5). In this case, we 
used importance sampling with the following Jastrow-like guiding state [9] 



,[0,t)N 



{n\g) = exp 



i,jeA 



Uj 

2 



(65) 



where = {xi,yi), a is a real positive parameter, and the long range potential v is 
defined as 



v[r 



Al ^ 



JqxX+iqyy 



1 -|- (cos Qx + COS qy)/2 
1 — (cos Qx + COS qy)/2 



(66) 
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Figure 7. Local energy per site lE{Mg,no'Hg{nt))/E{A4gjTig)]/ \M versus the 
imaginary time t for a 6 x 6 Heisenberg system with — ^^'^ J = 1. The 

Monte Carlo simulation (dots with error bars) was done by using importance sampling 
with the guiding function (67) and statistical parameters M — 2^^, N — 2^, and 
i? = 20. 



the sum over and qy being extended over the Brillouin zone 27r/L, 47r/L, . . . , 27r, with 
excluded. From Eq. (65) we have 
{n®lk®li\g) 



{n\g) 



exp 



a 



( - i ) [{uk © 1 - nfc) v{r, - Tk) 



jgA, i^k, I 

+ {ni®l- ni)v{ri - ri)] 



(67) 



We assumed a = 1.2 as suggested in Ref. [9]. 

In Fig. 7 we show the local energy per site [E{M.f,noTi.g{nt))/E{Mf,'^)]/ |A| as a 
function of the imaginary time t for a 6 x 6 Heisenberg system having X]l=i ~ 0- 
The amplitude of the error bars shown in Fig. 7 is considerably reduced with respect 
to the value that one would obtain without using importance sampling. We also stress 
that the dynamics shown in 7 is relative to the Hamiltonian Hg modified by the chosen 
guiding function \g). 

In Fig. 8 we show the local staggered magnetization {3[E{A4g%^ S^^{nt)Aif,'^) 
/E{Ai^g'^^'^)]/\A\Y^'^ evaluated in Heisenberg systems of different size as a function 
of the imaginary time t' and for a large value of the other imaginary time t. Here 
S^^{n) = {n\S^^\n) is the quantum expectation in the Fock state n of the diagonal 
operator 



5f 



1 

lAI 



{61 



As noticed in the case of Fig. 7, also the dynamics of the local staggered magnetization 
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Figure 8. Local staggered magnetization {3[E{A4g.nt 
/|A|}^/^, defined in terms of the operator (68), versus the imaginary time t' for 
different-size Heisenberg systems with ^ and J = 1. The Monte Carlo 

simulations (dots with error bars) were done by using importance sampling with the 
guiding function (67) with a = 1.2 and parameters t = 5, M = 2^^, N = 22, and 
i? = 10 in the 4 X 4 case, t = 10, M = 2^^, N = 22, and i? = 20 in the 6 x 6 case, 
t = 20, M = 2", N = 22, and i? = 40 in the 8 x 8 case. 



shown in Fig. 8 is relative to the Hamiltonian Hg modified by the guiding function 
(67). The asymptotic values of the local staggered magnetization reached for t' large 
are in agreement with those obtained with different Monte Carlo algorithms [9, 16]. The 
statistical errors shown in Fig. 8 can be reduced by a factor about 10 by exploiting the 
covariance between the local estimators for S^^^ and Eq, as explained in [16]. 

Finally, in Fig. 9 we provide an example of how the local expectation values shown 
in Fig. 8 depend on the parameter a of the guiding function (67). For different values 
of a the local expectations have different evolutions determined by Hg(a)^ however, as 
stated in Section 5.3 for a general guiding function, they all converge to the quantum 
expectation of S^^ in the ground state of H. In agreement with Ref. [9], the value 
a = 1.2 is close to the optimal choice, which provides smallest fluctuations and minimal 
evolution with respect to the asymptotic values. 

In Figs. 1 and 2 we have shown that the imaginary- and real-time evolution of 
the expectation of the basic functional coincides with that of the corresponding 

quantum matrix element Xln(^l^~^*l^o)- Of course, a similar behavior is general. 
Even if not shown explicitly, in all the considered examples the evolution of the relevant 
time-dependent probabilistic expectations coincides with that of the corresponding time- 
dependent quantum correlations functions. 



Exact Monte Carlo time dynamics in many-body lattice quantum systems 



+ s 



0.55 



23 



o 'v 



CO 




0.45 



Figure 9. Local staggered magnetization versus the imaginary time t' for the 4x4 
Heisenberg system of Fig. 8. The three curves were obtained with different values of the 
parameter a defining the importance sampling function (67). The other parameters 
were t = 5, M = 2^^, N = 22, and R = 20. 



8. Conclusions 

We have exploited an exact probabilistic representation of the quantum dynamics in 
a lattice to derive a Monte Carlo algorithm, named EPRMC, for which standard 
fluctuation control techniques like reconfigurations and importance sampling have been 
adapted and rigorously proved. This exact representation holds for both imaginary and 
real times, even if in the latter case only a partial fluctuation control is possible so that 
reliable statistical simulations are limited to short times. 

Monte Carlo algorithms, like GFMC or GFMCP, provide similar representations of 
the evolution operator, which are affected, however, by a systematic error e controlled 
by the number of iterations performed. With respect to these approximated methods, 
EPRMC gives an efficiency gain proportional to the accuracy 6~^. 
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Appendix 

In this Appendix, we calculate the relative efficiencies of GFMC and EPRMC methods. 
Both the methods have the aim to sample the operator e~^* for t large. We can write 

U{t) = e~"' ~ 6"^°*, for t > t, (69) 
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where t is the characteristic time to filter out the excited states Ei, E2, . . . with respect 
to the ground state £"0, 

El — Eq ^ ^ 

As explained in the introduction, GFMC samples the operator (1 — Ht/N)'^ whereas 
EPRMC samples directly the operator e"^*. Since limAr^oo(l + x/iV)^ = e^, GFMC 
— > EPRMC as the number of iterations in the GFMC method grows. However, for 
a finite value of A^, GFMC remains affected by a systematic error. We are interested in 
evaluating the critical value of A^ above which this error becomes smaller than a given 
value. Let us consider the relative difference 

fN{x) = ^'^^^ ■ (71) 



By using 



log(l + y) = J2 ^^2) 



k=l 



.2 . ^3 



Eq. (71) becomes 

/^(x) = (l-e-^+5^--). (73) 

For concreteness, let us put x = —E^t in Eq. (73). If we require that the relative error 
is fN{—Eot) = e 1, then we must have A^ > Nt{e), where 

Ntie) = (74) 

In conclusion, Nt{e) is the number of steps needed in GFMC to sample the operator 
Q-Ht ^ large with a relative error equal to e. On the other hand, the number of 
steps needed in EPRMC to sample e~^* for t large is given by the average number of 
jumps that, when an optimal reconfiguration scheme is chosen as discussed in Section 
7, coincides with the number of reconfigurations Rt 

Rt = {A)r]t ~ Eft, (75) 

where {A) is the average number of active links and E^l^^ is the ground-state energy in 
the non-interacting case. Therefore, the relative efficiency of EPRMC with respect to 
GFMC is given by the ratio 

Nt{e) Elt 



Rt 2Ef\' 



(76) 



We see that the superiority of EPRMC grows by increasing the time t or increasing the 
accuracy required in GFMC. In particular for t = t we have 



Ntie) E, 



- f77) 

Rt 2Ef\E,-E,)e 

It is clear that if, instead of GFMC, we consider the GFMCP method, the efficiency 
ratio (76) changes. In fact, any step in GFMCP, on the average, amounts to (n^) 
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elementary GFMC steps, where roughly {ug) = {A)rit [4]. Thus, in GFMCP the number 
of steps needed to sample the operator e~^* for t large with a relative error e is reduced to 
Nt{e) = Elt/{2Ef\) so that the relative efficiency of EPRMC with respect to GFMCP 
is given by the ratio 



required in GFMCP. 
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